Project

Mapping of Ki67 interactions with the genome and comparison with lamina interactions.

Introduction

Various analyses of RPE cells synchronized in mitosis andwith time points during the cell cycle afterwards.

Method

NA

Set-up

Set the parameters and list the data.

# Load dependencies - without warnings / messages
library(tidyverse)
library(GenomicRanges)
library(rtracklayer)
library(ggplot2)
library(RColorBrewer)
library(GGally)
library(corrr)
library(caTools)
library(ggbeeswarm)

# Prepare output 
output_dir <- "ts210519_ki67_rpe_cell_cycle"
dir.create(output_dir, showWarnings = FALSE)


# Load input
chromosomes <- c(paste0("chr", 1:22), "chrX")


input_dir <- "ts210413_data_gathering"

bin_size <- readRDS(file.path(input_dir, "bin_size.rds"))
centromeres <- readRDS(file.path(input_dir, "centromeres.rds"))

colors_set1 <- readRDS(file.path(input_dir, "colors_set1.rds"))
colors_set2 <- readRDS(file.path(input_dir, "colors_set2.rds"))
colors_set3 <- readRDS(file.path(input_dir, "colors_set3.rds"))

tib_padamid_replicates <- readRDS(file.path(input_dir, 
                                            "tib_padamid_replicates.rds"))
tib_padamid_combined <- readRDS(file.path(input_dir, 
                                          "tib_padamid_combined.rds"))

gr_padamid_replicates <- readRDS(file.path(input_dir, 
                                           "gr_padamid_replicates.rds"))
gr_padamid_combined <- readRDS(file.path(input_dir, 
                                         "gr_padamid_combined.rds"))

tib_hmm_replicates <- readRDS(file.path(input_dir, "tib_hmm_replicates.rds"))
tib_hmm_combined <- readRDS(file.path(input_dir, "tib_hmm_combined.rds"))

padamid_metadata_replicates <- readRDS(file.path(input_dir, 
                                                 "padamid_metadata_replicates.rds"))
padamid_metadata_combined <- readRDS(file.path(input_dir, 
                                               "padamid_metadata_combined.rds"))



# Prepare seqnames
chrom_sizes <- tibble(seqnames = seqlevels(gr_padamid_combined),
                      length = seqlengths(gr_padamid_combined)) %>%
  arrange(-length)


# Scale pA-DamID scores?
tib_padamid_combined <- tib_padamid_combined %>%
  mutate_at(4:ncol(.), function(x) scale(x)[, 1])
library(knitr)
opts_chunk$set(fig.width = 10, fig.height = 4, cache = T,
               dev=c('png', 'pdf'), fig.path = file.path(output_dir, "figures/")) 
pdf.options(useDingbats = FALSE)
# From Fede:
# ggpairs custom functions
corColor <- function(data, mapping, color = I("black"), sizeRange = c(1, 3), ...) {

  x   <- eval_data_col(data, mapping$x)
  y   <- eval_data_col(data, mapping$y)
  r   <- cor(x, y)
  rt  <- format(r, digits = 3)
  tt  <- as.character(rt)
  cex <- max(sizeRange)

  # helper function to calculate a useable size
  percent_of_range <- function(percent, range) {
    percent * diff(range) + min(range, na.rm = TRUE)
  }

  # plot correlation coefficient
  p <- ggally_text(label = tt, mapping = aes(), xP = 0.5, yP = 0.5,
                   # size = I(percent_of_range(cex * abs(r), sizeRange)), 
                   size = 6, 
                   color = color, ...) +
    theme(panel.grid.minor=element_blank(),
          panel.grid.major=element_blank())

  corColors <- RColorBrewer::brewer.pal(n = 7, name = "RdYlBu")[2:6]

  if (r <= boundaries[1]) {
    corCol <- corColors[1]
  } else if (r <= boundaries[2]) {
    corCol <- corColors[2]
  } else if (r < boundaries[3]) {
    corCol <- corColors[3]
  } else if (r < boundaries[4]) {
    corCol <- corColors[4]
  } else {
    corCol <- corColors[5]
  }

  p <- p +
    theme(panel.background = element_rect(fill = corCol))

  return(p)
}

customScatter <- function (data, mapping) 
{
    p <- ggplot(data = data, mapping) + 
      geom_bin2d(bins = 100) +
      geom_smooth(method = "lm", se = T, col = "red") +
      scale_fill_gradient(low = "lightgrey", high = "black", name = "Count") +
      theme_bw()
    
    p 
}

PlotScatter <- function(tib, n1, n2, color_by = NULL, identity = F,
                        xlimits = NULL, facet_seqnames = F,
                        smooth_line = F, smooth_seqnames = F) {
  # Get tibble
  tib <- tib %>%
    dplyr::select(seqnames, matches(n1), matches(n2)) %>%
    rename_all(~ c("seqnames", "n1", "n2")) %>%
    mutate(seqnames = factor(seqnames, 
                             levels = seqlevels(gr_padamid_combined)))
  
  # Prepare color
  if (! is.null(color_by)) {
    tib <- tib %>%
      add_column(color = color_by) %>%
      drop_na()
    alpha = 1
    limits_color <- quantile(tib$color, c(0.001, 0.999), na.rm = T)
    tib$color[tib$color < limits_color[1]] <- limits_color[1]
    tib$color[tib$color > limits_color[2]] <- limits_color[2]
  } else {
    tib <- tib %>% drop_na()
    tib$color = "1"
    alpha = 0.02
  }
  
  # Plot
  if (is.null(xlimits)) {
    xlimits <- quantile(tib$n1, c(0.001, 0.999), na.rm = T) * 1.4 
  }
  ylimits <- quantile(tib$n2, c(0.001, 0.999), na.rm = T) * 1.4
  
  plt <- tib %>%
    arrange(sample(1:nrow(.), size = nrow(.), replace = F)) %>%
    ggplot(aes(x = n1, y = n2, color = color)) +
    geom_point(size = 0.5, alpha = alpha) +
    geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
    geom_vline(xintercept = 0, col = "black", linetype = "dashed") +
    xlab(n1) +
    ylab(n2) +
    ggtitle(paste0("Spearman: ", 
                   round(cor(tib$n1, tib$n2, use = "complete.obs",
                             method = "spearman"), 2))) +
    coord_cartesian(xlim = xlimits, ylim = ylimits) +
    theme_bw() +
    theme(aspect.ratio = 1)
  
  # Prepare color
  if (! is.null(color_by)) {
    plt <- plt +
      scale_color_gradient2(low = "blue", mid = "grey", high = "red",
                            midpoint = 0)
  } else {
    plt <- plt + 
      scale_color_manual(values = "black", guide = F)
  }
  if (identity) plt <- plt + geom_abline(slope = 1, intercept = 0, col = "red", linetype = "dashed")
  
  # Facets / smooth
  if (facet_seqnames) plt <- plt + facet_wrap(~ seqnames)
  if (smooth_line) plt <- plt + geom_smooth(method = "loess", span = 0.7, se = F, col = "black")
  if (smooth_seqnames) plt <- plt + geom_smooth(aes(group = seqnames), 
                                                alpha = 0.2, col = "black",
                                                method = "loess", span = 0.7, se = F)
  
  plot(plt)
  
}

PlotScatterBinned <- function(tib, n1, n2, color_by = NULL, identity = F, 
                              n_min = 10, ylimits_col = c(-2.4, 2.4),
                              xlimits = c(0, 50), smooth_line = F,
                              limits_count = c(0, 400)) {
  # Get tibble
  tib_process <- tib %>%
    dplyr::select(seqnames, all_of(n1), all_of(n2)) %>%
    rename_all(~ c("seqnames", "n1", "n2"))
  
  if (! is.null(color_by)) {
    tib_process <- tib_process %>%
      add_column(color = color_by)
  }
  
  tib_process <- tib_process %>%
    drop_na() %>%
    # Filter xlimits before making the bins
    filter(n1 < (xlimits[2]+5))
  
  # Change color range
  if (! is.null(color_by)) {
    limits_color <- quantile(tib_process$color, c(0.001, 0.999), na.rm = T)
    tib_process$color[tib_process$color < limits_color[1]] <- limits_color[1]
    tib_process$color[tib_process$color > limits_color[2]] <- limits_color[2]
  }
  
  # Metrics
  n1_min = min(tib_process$n1) - 0.001
  n1_max = max(tib_process$n1) + 0.001
  n1_binsize <- (n1_max - n1_min) / 49
  
  n2_min = min(tib_process$n2) - 0.001
  n2_max = max(tib_process$n2) + 0.001
  n2_binsize <- (n2_max - n2_min) / 49
  
  tib_summary <- tib_process %>%
    mutate(n1_cut = cut(n1, 
                        seq(n1_min, n1_max, length.out = 50)),
           n2_cut = cut(n2, 
                        seq(n2_min, n2_max, length.out = 50))) %>%
    mutate(n1_bin = as.numeric(as.factor(n1_cut)),
           n2_bin = as.numeric(as.factor(n2_cut))) %>%
    mutate(n1_bin = n1_min - n1_binsize/2 + n1_bin * n1_binsize,
           n2_bin = n2_min - n2_binsize/2 + n2_bin * n2_binsize) %>%
    group_by(n1_bin, n2_bin)
  
  # Plot
  if (! is.null(color_by)) {
    tib_summary <- tib_summary %>%
    dplyr::summarise(n = n(),
                     mark = mean(color)) %>%
    ungroup() %>%
    filter(n >= n_min)
    
    plt <- tib_summary %>%
      ggplot(aes(x = n1_bin, y = n2_bin)) +
      geom_tile(aes(fill = mark)) +
      scale_fill_gradient2(low = "blue", mid = "grey", high = "red",
                           midpoint = 0, limits = ylimits_col, 
                           na.value = "green") +
      geom_hline(yintercept = 0, linetype = "dashed", col = "black")
  } else {
    tib_summary <- tib_summary %>%
    dplyr::summarise(n = n()) %>%
    ungroup() %>%
    filter(n >= n_min)
    
    plt <- tib_summary %>%
      ggplot(aes(x = n1_bin, y = n2_bin)) +
      geom_tile(aes(fill = n)) +
      scale_fill_gradient(low = "lightgrey", high = "black", name = "Count",
                          limits = limits_count, na.value = "green") +
      geom_hline(yintercept = 0, linetype = "dashed", col = "blue")
  }
  
  plt <- plt + 
    xlab(n1) +
    ylab(n2) +
    coord_cartesian(xlim = xlimits) +
    ggtitle(paste0("Pearson: ", 
                   round(cor(tib_process$n1, tib_process$n2, use = "complete.obs",
                             method = "pearson"), 2))) +
    theme_bw() +
    theme(aspect.ratio = 1)
  
  if (smooth_line) plt <- plt + geom_smooth(data = tib_process,
                                            aes(x = n1, y = n2),
                                            method = "loess", 
                                            span = 0.7, se = F, col = "red")
  
  
  plot(plt)
  
}


PlotDataTracksLight <- function(input_tib, exp_names, centromeres, 
                                color_groups, plot_chr = "chr1", 
                                plot_start = 1, plot_end = 500e6,
                                smooth = 1, color_list = NULL,
                                fix_yaxis = F) {
  
  # Exp names is a vector of sample names
  exp_search <- paste(exp_names, collapse = "|")
  
  # Get the scores for the samples
  tib_plot <- input_tib %>%
    dplyr::select(seqnames, start, end, 
                  matches(exp_search))
  
  if (smooth > 1) {
    tib_plot <- tib_plot %>%
      mutate_at(vars(contains("_")), 
                runmean, k = smooth)
  }
  
  # Filter for plotting window
  tib_plot <- tib_plot %>%
    filter(seqnames == plot_chr,
           start >= plot_start,
           end <= plot_end)
  
  # Gather
  tib_gather <- tib_plot %>%
    gather(key, value, -seqnames, -start, -end) %>%
    mutate(key = factor(key, levels = exp_names),
           fill_column = color_groups[match(key, names(color_groups))],
           fill_column = factor(fill_column, levels = unique(color_groups)))
  
  
  # Centromeres
  centromeres.plt <- as_tibble(centromeres) %>%
    filter(seqnames == plot_chr) %>%
    gather(key_centromeres, value, start, end)
  
  # Plot
  ylimits <- quantile(tib_gather$value, c(0.001, 0.999), na.rm = T)
  fill_column <- NULL
  
  plt <- tib_gather %>% 
    ggplot(aes(fill = fill_column))
  
  
  # Set all counts tracks to the same limits
  plt <- plt + 
    geom_rect(aes(xmin = start / 1e6, xmax = end / 1e6, 
                  ymin = 0, ymax = value)) +
    geom_line(data = centromeres.plt, 
              aes(x = value / 1e6, y = 0),
              color = "black", size = 2) +
    geom_hline(yintercept = 0, size = 0.5) +
    facet_grid(key ~ ., scales = "free_y") +
    xlab(paste0(plot_chr, " (Mb)")) +
    ylab("Score") +
    scale_x_continuous(expand = c(0, 0)) + 
    scale_y_continuous(expand = c(0.025, 0.025)) +
    theme_classic()
  
  if (! is.null(color_list)) {
    colors <- levels(tib_gather$fill_column)

    color_list <- color_list[1:length(colors)]
    names(color_list) <- colors
    #colors <- recode(colors, !!!color_list)
    
    plt <- plt +
      scale_fill_manual(values = color_list, guide = F)
  } else {
    plt <- plt +
      scale_fill_brewer(palette = "Set1", guide = F)
  }
  
  if (fix_yaxis) {
    plt <- plt + 
      coord_cartesian(xlim = c(max(c(plot_start,
                                   min(tib_gather$start))) / 1e6,
                             min(c(plot_end,
                                   max(tib_gather$end))) / 1e6),
                    ylim = ylimits)
  } else {
    plt <- plt + 
      coord_cartesian(xlim = c(max(c(plot_start,
                                   min(tib_gather$start))) / 1e6,
                             min(c(plot_end,
                                   max(tib_gather$end))) / 1e6))
  }
  
  plot(plt)
  
}

quantiles <- function(x) {
  # Use quantiles as boxplot boundaries
  r <- quantile(x, probs = c(0.05, 0.25, 0.5, 0.75, 0.95))
  names(r) <- c("ymin", "lower", "middle", "upper", "ymax")
  r
}

1. Correlations between time points

I will perform these analyses only on the combined data, as the data is only shallowly sequenced.

# Gather metadata
padamid_metadata_cellcycle <- padamid_metadata_combined %>%
  filter(cell == "RPE" &
           experiment %in% c("wildtype", "cell_cycle"),
         condition != "t_21h")
         #! target %in% c("LMNB1", "H3K27me3", "H3K9me3"))

# Use GGally to make correlation plots
boundaries <- seq(from = 0.1, to = 0.7, length.out = 4)

# Combined
tib <- tib_padamid_combined %>%
  drop_na() %>%
  dplyr::select(padamid_metadata_cellcycle$sample)

plt <- ggpairs(tib %>%
                 dplyr::select(contains("Ki67")),
               upper = list(continuous = corColor),
               lower = list(continuous = customScatter),
               diag = list(continuous = function(data, mapping, ...) {
                   ggally_densityDiag(data = data, mapping = mapping, alpha = 0.3, fill = "red") +
                   theme_bw()})) +
  ggtitle("Correlating all data") +
  xlab("") +
  ylab("")

print(plt)

# Correlation plot only
metadata <- padamid_metadata_cellcycle

tib_cor <- tib %>%
  dplyr::select(contains("Ki67")) %>%
  correlate(method = "pearson") %>%
  gather(key, value, -term) %>%
  mutate(value = replace_na(value, 1)) %>%
  mutate(n1 = str_remove(term, "_Ki67"),
         n2 = str_remove(key, "_Ki67"),
         match1 = match(term, metadata$sample),
         match2 = match(key, metadata$sample),
         experiment1 = metadata$experiment[match1],
         experiment2 = metadata$experiment[match2],
         condition1 = metadata$condition[match1],
         condition2 = metadata$condition[match2]) %>%
  drop_na()

# Plot
plt <- tib_cor %>%
  ggplot(aes(x = condition1, y = condition2, fill = value)) +
  geom_tile() +
  xlab("") + ylab("") +
  scale_fill_gradientn(colors = colorRampPalette(rev(brewer.pal(n = 7, 
                                                                name = "RdYlBu")))(100),
                       limits = c(min(0, tib_cor$value), 1),
                       name = "Pearson correlation") +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, hjust = 1))
plot(plt)

As expected, t = 0h does not correlate well with the other time points. Also, correlations get lower if time points are further apart.

I also want to make these plots for individual chromosomes. Scatterplots are too much, but simple pearson correlations should do the trick.

# Calculate pearson correlations per chromosome
CorrelationsPerChromosome <- function(input_tib, metadata) {
  
  # Gather data and calculate pearson correlations
  tib <- input_tib %>%
    drop_na() %>%
    dplyr::select(seqnames, metadata$sample)
  
  # Calculate mean score per chromosome
  tib_summary <- tib %>%
    gather(key, value, -seqnames) %>%
    mutate(match = match(key, metadata$sample),
           experiment = metadata$experiment[match],
           condition = metadata$condition[match],
           seqnames = factor(seqnames, chromosomes)) %>%
    mutate(size = seqlengths(gr_padamid_combined)[match(seqnames, 
                                                        seqlevels(gr_padamid_combined))],
           size = size / 1e6) %>%
    drop_na() 
  
  # Boxplot of all scores per chromosome
  plt <- tib_summary %>%
    arrange(-size) %>% 
    mutate(seqnames = factor(seqnames, levels = unique(seqnames))) %>%
    ggplot(aes(x = seqnames, y = value, fill = condition)) +
    #geom_boxplot(outlier.shape = NA) +
    stat_summary(fun.data = quantiles, geom = "boxplot") +
    geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
    facet_grid(. ~ condition) +
    xlab("") +
    ylab("Ki67 score") +
    coord_cartesian(ylim = c(-2.5, 3.5)) +
    scale_fill_manual(values = colors_set2, guide = F) +
    theme_bw() +
    theme(aspect.ratio = 1,
          axis.text.x = element_text(angle = 90, hjust = 1))
  plot(plt)
  
  plt <- tib_summary %>%
    arrange(-size) %>% 
    mutate(seqnames = factor(seqnames, levels = unique(seqnames))) %>%
    filter(condition %in% c("wt", "t_1h")) %>%
    ggplot(aes(x = seqnames, y = value, fill = condition)) +
    #geom_boxplot(outlier.shape = NA) +
    stat_summary(fun.data = quantiles, geom = "boxplot",
                 position = "dodge") +
    geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
    xlab("") +
    ylab("Ki67 score") +
    coord_cartesian(ylim = c(-2.5, 3.5)) +
    scale_fill_manual(values = c("grey30", "#E41A1C")) +
    theme_bw() +
    theme(aspect.ratio = 1,
          axis.text.x = element_text(angle = 90, hjust = 1))
  plot(plt)
  
  
  # Boxplot of all scores per chromosome - by condition
  plt <- tib_summary %>%
    arrange(-size) %>% 
    filter(! condition %in% c("wt", "t_0h")) %>%
    mutate(seqnames = factor(seqnames, levels = unique(seqnames))) %>%
    ggplot(aes(x = condition, y = value, fill = condition)) +
    #geom_boxplot(outlier.shape = NA) +
    stat_summary(fun.data = quantiles, geom = "boxplot") +
    geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
    facet_wrap( ~ seqnames) +
    xlab("") +
    ylab("Ki67 score") +
    coord_cartesian(ylim = c(-1.5, 1.5)) +
    scale_fill_manual(values = colors_set2, guide = F) +
    theme_bw() +
    theme(aspect.ratio = 1,
          axis.text.x = element_text(angle = 90, hjust = 1))
  plot(plt)
  
  
  # Single point (mean score) per chromosome - without wildtype and t=0h
  plt <- tib_summary %>%
    group_by(experiment, condition, size, seqnames) %>%
    summarise(mean = mean(value)) %>%
    filter(! condition %in% c("wt", "t_0h")) %>%
    ungroup() %>% 
    arrange(-size) %>% 
    mutate(seqnames = factor(seqnames, levels = unique(seqnames))) %>%
    ggplot(aes(x = seqnames, y = mean, col = condition)) +
    geom_point(size = 2) +
    geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
    xlab("") +
    ylab("Mean Ki67 score") +
    scale_color_brewer(palette = "Spectral") +
    theme_bw() +
    theme(aspect.ratio = 1,
          axis.text.x = element_text(angle = 90, hjust = 1))
  plot(plt)
  
  # Same, but difference plot
  plt <- tib_summary %>%
    group_by(experiment, condition, size, seqnames) %>%
    summarise(mean = mean(value)) %>%
    filter(! condition %in% c("wt", "t_0h")) %>%
    ungroup() %>% 
    spread(condition, mean) %>%
    mutate(t_3h = t_3h - t_1h,
           t_6h = t_6h - t_1h,
           t_10h = t_10h - t_1h) %>%
    dplyr::select(-t_1h) %>%
    gather(condition, mean, contains("t_")) %>%
    arrange(-size) %>% 
    mutate(seqnames = factor(seqnames, levels = unique(seqnames)),
           condition = factor(condition, levels = levels(metadata$condition))) %>%
    ggplot(aes(x = seqnames, y = mean, col = condition)) +
    geom_point(size = 2) +
    geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
    xlab("") +
    ylab("Difference in Ki67 with t=1h") +
    scale_color_manual(values = brewer.pal("Spectral", n = 4)[2:4]) +
    theme_bw() +
    theme(aspect.ratio = 1,
          axis.text.x = element_text(angle = 90, hjust = 1))
  plot(plt)
  
  
  # Can I show that interactions are increased on large chromosomes compared to unsync cells?
  plt <- tib %>%
    mutate(t_1h = RPE_1h_Ki67 - RPE_wt_Ki67,
           t_3h = RPE_3h_Ki67 - RPE_wt_Ki67,
           t_6h = RPE_6h_Ki67 - RPE_wt_Ki67,
           t_10h = RPE_10h_Ki67 - RPE_wt_Ki67) %>%
    dplyr::select(seqnames, starts_with("t_")) %>%
    gather(key, value, -seqnames) %>%
    mutate(seqnames = factor(seqnames, chrom_sizes$seqnames),
           key = factor(key, levels = c("t_1h", "t_3h",
                                        "t_6h", "t_10h"))) %>%
    ggplot(aes(x = seqnames, y = value)) +
    #geom_boxplot(outlier.shape = NA) +
    stat_summary(fun.data = quantiles, geom = "boxplot") +
    facet_grid(. ~ key) +
    theme_bw() +
    theme(aspect.ratio = 1)
  #plot(plt)
  
  # By chromosome
  plt <- tib %>%
    mutate(t_1h = RPE_1h_Ki67 - RPE_wt_Ki67,
           t_3h = RPE_3h_Ki67 - RPE_wt_Ki67,
           t_6h = RPE_6h_Ki67 - RPE_wt_Ki67,
           t_10h = RPE_10h_Ki67 - RPE_wt_Ki67) %>%
    dplyr::select(seqnames, starts_with("t_")) %>%
    gather(key, value, -seqnames) %>%
    mutate(seqnames = factor(seqnames, chrom_sizes$seqnames),
           key = factor(key, levels = c("t_1h", "t_3h",
                                        "t_6h", "t_10h"))) %>%
    filter(key == "t_1h") %>%
    ggplot(aes(x = seqnames, y = value, fill = seqnames)) +
    #geom_boxplot(outlier.shape = NA) +
    stat_summary(fun.data = quantiles, geom = "boxplot") +
    geom_hline(yintercept = 0, col = "grey30", linetype = "dashed") +
    coord_cartesian(ylim = c(-2.3, 1.5)) +
    xlab("Chromosome (by size)") +
    ylab("Difference in Ki67 with unsync") +
    scale_fill_discrete(guide = F) +
    theme_bw() +
    theme(aspect.ratio = 1,
          axis.text.x = element_text(angle = 90, hjust = 1))
  plot(plt)
    
  
  # Fraction of the genome "enriched"
  plt <- tib %>%
    drop_na() %>%
    gather(key, value, -seqnames) %>%
    group_by(key, seqnames) %>%
    dplyr::summarise(fraction = mean(value > 0)) %>%
    ungroup() %>%
    mutate(seqnames = factor(seqnames, chrom_sizes$seqnames),
           key = factor(key, names(tib))) %>%
    filter(key %in% c("RPE_wt_Ki67", "RPE_1h_Ki67")) %>%
    ggplot(aes(x = seqnames, y = fraction, fill = key)) +
    geom_bar(stat = "identity", position = "dodge") +
    scale_fill_manual(values = c("grey30", "#E41A1C")) +
    theme_bw() +
    theme(aspect.ratio = 1,
          axis.text.x = element_text(angle = 90, hjust = 1))
  #plot(plt)
  
  
  # Mean score vs chromosome size
  plt <- tib_summary %>%
    group_by(experiment, condition, size, seqnames) %>%
    summarise(mean = mean(value)) %>%
    ggplot(aes(x = size, y = mean, col = condition)) +
    geom_point() +
    geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
    facet_grid(. ~ condition) +
    xlab("Chromosome size (Mb)") +
    ylab("Mean Ki67 score") +
    coord_cartesian(xlim = c(0, 250)) +
    scale_color_manual(values = colors_set2, guide = F) +
    theme_bw() +
    theme(aspect.ratio = 1,
          axis.text.x = element_text(angle = 90, hjust = 1))
  plot(plt)
  
  tib_cor <- tib %>%
    group_by(seqnames) %>%
    nest() %>%
    mutate(cordf = map(data, correlate)) %>%
    unnest(cordf) %>%
    dplyr::select(-data) %>%
    gather(key, value, -seqnames, -term) %>%
    drop_na() %>%
    mutate(n1 = str_remove(term, "_Ki67"),
           n2 = str_remove(key, "_Ki67"),
           match1 = match(term, metadata$sample),
           match2 = match(key, metadata$sample),
           experiment1 = metadata$experiment[match1],
           experiment2 = metadata$experiment[match2],
           condition1 = metadata$condition[match1],
           condition2 = metadata$condition[match2],
           seqnames = factor(seqnames, chromosomes)) %>%
    mutate(size = seqlengths(gr_padamid_replicates)[match(seqnames, 
                                                          seqlevels(gr_padamid_replicates))],
           size = size / 1e6) %>%
    drop_na()
  
  # Plot
  plt <- tib_cor %>%
    arrange(-size) %>% 
    mutate(seqnames = factor(seqnames, levels = unique(seqnames))) %>%
    ggplot(aes(x = seqnames, y = value, col = condition2)) +
    geom_point() +
    geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
    facet_grid(. ~ condition1) +
    xlab("") +
    ylab("Pearson correlation") +
    coord_cartesian(ylim = c(0, 1)) +
    scale_color_manual(values = colors_set2) +
    theme_bw() +
    theme(aspect.ratio = 1,
          axis.text.x = element_text(angle = 90, hjust = 1))
  plot(plt)
  
  plt <- tib_cor %>%
    ggplot(aes(x = size, y = value, col = condition2)) +
    geom_point() +
    geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
    facet_grid(. ~ condition1) +
    xlab("Chromosome size (Mb)") +
    ylab("Pearson correlation") +
    coord_cartesian(xlim = c(0, 250),
                    ylim = c(0, 1)) +
    scale_color_manual(values = colors_set3) +
    theme_bw() +
    theme(aspect.ratio = 1,
          axis.text.x = element_text(angle = 90, hjust = 1))
  plot(plt)
  
}

# Correlations per chromosome
CorrelationsPerChromosome(tib_padamid_combined, 
                          padamid_metadata_cellcycle %>%
                            filter(target == "Ki67"))

2. Enrichment at centromeres and telomeres

These analyses are based on simple counts. I want to do a similar analysis based on the normalized scores themselves.

# Get similar chromosomes for the two objects and filter samples
tib_padamid_cellcycle <- tib_padamid_combined %>%
  dplyr::select(seqnames, start, end, padamid_metadata_cellcycle$sample) %>%
  filter(seqnames %in% chromosomes)
centromeres <- centromeres[seqnames(centromeres) %in% chromosomes]

# Also, prepare telomeres
telomeres <- c(GRanges(seqnames = chromosomes,
                       ranges = IRanges(start = 1, end = 1)),
               GRanges(seqnames = chromosomes,
                       ranges = IRanges(start = seqlengths(gr_padamid_combined)[chromosomes], 
                                        end = seqlengths(gr_padamid_combined)[chromosomes])))
seqinfo(telomeres) <- seqinfo(gr_padamid_combined)
telomeres <- sort(telomeres)

# Add distance to centromeres and telomeres to the GRanges object - in Mb
dis <- distanceToNearest(as(tib_padamid_cellcycle, "GRanges"), centromeres)
tib_padamid_cellcycle$distance_to_centromere <- mcols(dis)$distance / 1e6

dis <- distanceToNearest(as(tib_padamid_cellcycle, "GRanges"), telomeres)
tib_padamid_cellcycle$distance_to_telomere <- mcols(dis)$distance / 1e6

# Process
tib <- tib_padamid_cellcycle %>%
  dplyr::select(-start, -end) %>%
  mutate(rdna = as.logical(seqnames %in% paste0("chr", c(13, 14, 15, 21, 22)))) %>%
  gather(key, value, -contains("distance"), -rdna, -seqnames) %>%
  separate(key, c("cell", "condition", "target"), remove = F) %>%
  filter(target == "Ki67") %>%
  mutate(match = match(key, padamid_metadata_cellcycle$sample),
         condition = padamid_metadata_cellcycle$condition[match],
         experiment = padamid_metadata_cellcycle$experiment[match],
         dis_group_centromere = as.numeric(cut(distance_to_centromere, 
                                               breaks = seq(-1, max(distance_to_centromere) + 5, 
                                                            by = 1))) - 1,
         dis_group_telomere = as.numeric(cut(distance_to_telomere, 
                                             breaks = seq(-1, max(distance_to_telomere) + 5, 
                                                          by = 1))) - 2,
         dis_group_telomere = ifelse(dis_group_telomere < 0, 0, dis_group_telomere)) %>%
  dplyr::select(-contains("distance")) %>%
  gather(class, class_value, contains("dis_group")) %>%
  mutate(class = str_remove(class, "dis_group_"))


# Plot as boxplots
# First, calculate boxplots to have more control over ylimits
tib_summary <- tib %>%
  group_by(experiment, condition, class, class_value) %>%
  drop_na() %>%
  summarise(ymin = quantile(value, 0.05), 
            lower = quantile(value, 0.25), 
            middle = quantile(value, 0.5), 
            upper = quantile(value, 0.75), 
            ymax = quantile(value, 0.95))

tib_summary %>%
  filter(class_value <= 30) %>%
  ggplot(aes(x = class_value, ymin = ymin, lower = lower, middle = middle, 
             upper = upper, ymax = ymax, group = class_value, y = middle)) +
    geom_boxplot(stat="identity", outlier.shape = NA, fill = "lightgrey", col = "black") +
    geom_line(aes(y = middle, group = NULL), col = "red") +
    geom_hline(yintercept = 0, col = "blue", linetype = "dashed") +
    facet_grid(class ~ condition, scales = "free_y") + 
    xlab("Distance to feature (Mb)") +
    ylab("pA-DamID (log2)") +
    #coord_cartesian(xlim = c(0, 30)) +
    theme_bw() +
    theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, hjust = 1))

# Is the enrichment specific for rDNA chromosomes, or is chromosome size more 
# important?
tib_summary <- tib %>%
  filter(class_value < 2) %>%
  mutate(condition = droplevels(condition)) %>%
  group_by(experiment, condition, seqnames, class) %>%
  drop_na() %>%
  summarise(mean = mean(value)) %>%
  ungroup() %>%
  mutate(seqnames = factor(seqnames, chromosomes),
         rdna = seqnames %in% paste0("chr", c(13, 14, 15, 21, 22))) %>%
  drop_na() %>%
  mutate(size = seqlengths(gr_padamid_replicates)[match(seqnames, 
                                                        seqlevels(gr_padamid_replicates))],
         size = size / 1e6)

tib_summary %>%
  arrange(-size) %>% 
  mutate(seqnames = factor(seqnames, levels = unique(seqnames))) %>%
  ggplot(aes(x = seqnames, y = mean, col = class, shape = rdna)) +
  geom_point(size = 2) +
  geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
  facet_grid(. ~ condition) +
  xlab("") +
  ylab("Ki67 score near centromeres") +
  scale_color_manual(values = colors_set3) +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

# Add chromosome size and make a scatter plot
tib_summary %>%
  ggplot(aes(x = size, y = mean, col = class, shape = rdna)) +
  geom_point(size = 2) +
  geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
  xlab("Chromosome size (Mb)") +
  ylab("Ki67 score near centromeres (<2Mb)") +
  scale_color_manual(values = colors_set3) +
  facet_grid(. ~ condition) +
  coord_cartesian(xlim = c(0, 250)) +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, hjust = 1))

This mostly shows the expected result: t=0h is enriched near telomeres while the other time points are more enriched at/near centromeres.

3. Overlap in HMM calls

To determine the best strategy for the analysis, I need to consider a few options. I previously used HMM calls to do this. Can I also do this for the Ki67 data, or is this of too low quality? Let’s find out.

GenomeOverlap <- function(tib_hmm, metadata) {
  # Genome overlap is defined as the number of enriched domains (AD) over 
  # not-enriched (iAD) domains. NA bins are discarded
  
  tib <- tib_hmm %>%
    drop_na() %>%
    dplyr::select(metadata$sample) %>%   
    mutate_all(function(x) (x == "AD") + 0) %>%
    #rename_all(function(x) str_remove(x, "_Ki67")) %>%
    summarise_all(.funs = mean) %>%
    gather(key, value) %>%
    mutate(match = match(key, metadata$sample),
         condition = metadata$condition[match],
         experiment = metadata$experiment[match])
  
  plt <- tib %>%
    ggplot(aes(x = condition, y = value, fill = experiment)) +
    geom_bar(stat = "identity") +
    geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
    xlab("") +
    ylab("Overlap genome") +
    scale_fill_manual(values = colors_set3, guide = F) +
    theme_bw() +
    theme(aspect.ratio = 1,
          axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
  
  plot(plt)
  
}


PlotJaccard <- function(tib_hmm, metadata) {
  # Function to create heatmap of Jaccard indices
  
  tib <- tib_hmm %>%
    drop_na() %>%
    dplyr::select(metadata$sample) %>%   
    mutate_all(function(x) (x == "AD") + 0)
  
  
  # For this, use a correlation heatmap with heatmap 
  tib_dist <- 1 - as.matrix(stats::dist(t(as.matrix(tib)), method = "binary"))
  tib_dist <- as_tibble(tib_dist) %>%
    gather(key, value) %>%
    mutate(match = match(key, metadata$sample),
           condition = metadata$condition[match],
           experiment = metadata$experiment[match],
           target = rep(unique(condition), times = length(unique(condition))))
  
  plt <- tib_dist %>%
    ggplot(aes(x = condition, y = target, fill = value)) +
    geom_tile() +
    xlab("") + ylab("") +
    scale_fill_gradientn(colors = colorRampPalette(rev(brewer.pal(n = 7, 
                                                                  name = "RdYlBu")))(100),
                         limits = c(0, 1),
                         name = "Jaccard index") +
    theme_bw() +
    theme(aspect.ratio = 1,
          axis.text.x = element_text(angle = 90, hjust = 1))
  plot(plt)
  
}

# Plot overlap genome
GenomeOverlap(tib_hmm_combined, 
              padamid_metadata_cellcycle %>%
                filter(target == "Ki67"))

# Make Jaccard heatmap
PlotJaccard(tib_hmm_combined, 
            padamid_metadata_cellcycle %>%
                filter(target == "Ki67"))

Maybe I can use this for the analyses when the data is sequenced deeper. Right now, the quality is insufficient to get stable HMM calls.

5. Mitosis pattern in more detail

tib <- tib_padamid_cellcycle

# Scatter plots
PlotScatter(tib, "RPE_0h_Ki67", "RPE_1h_Ki67", identity = T)

PlotScatter(tib, "RPE_0h_Ki67", "RPE_wt_H3K27me3", identity = T)

PlotScatter(tib, "RPE_0h_Ki67", "RPE_wt_H3K9me3", identity = T)

PlotScatter(tib, "RPE_1h_Ki67", "RPE_wt_H3K27me3", identity = T)

PlotScatter(tib, "RPE_1h_Ki67", "RPE_wt_H3K9me3", identity = T)

# Combining telomeres distance and histone modifications 
PlotScatter(tib, "distance_to_centromere", "RPE_0h_Ki67", xlimits = c(0, 50))

PlotScatter(tib, "distance_to_centromere", "RPE_0h_Ki67", xlimits = c(0, 50), 
            color_by = tib$RPE_wt_H3K27me3)

PlotScatter(tib, "distance_to_centromere", "RPE_0h_Ki67", xlimits = c(0, 50), 
            color_by = tib$RPE_wt_H3K9me3)

PlotScatter(tib, "distance_to_telomere", "RPE_0h_Ki67", xlimits = c(0, 50))

PlotScatter(tib, "distance_to_telomere", "RPE_0h_Ki67", xlimits = c(0, 50), 
            color_by = tib$RPE_wt_H3K27me3)

PlotScatter(tib, "distance_to_telomere", "RPE_0h_Ki67", xlimits = c(0, 50), 
            color_by = tib$RPE_wt_H3K9me3)

PlotScatter(tib, "distance_to_telomere", "RPE_0h_Ki67", xlimits = c(0, 50), 
            smooth_line = T)
## `geom_smooth()` using formula 'y ~ x'

PlotScatter(tib, "distance_to_telomere", "RPE_0h_Ki67", xlimits = c(0, 50), 
            color_by = tib$RPE_wt_H3K27me3, smooth_line = T)
## `geom_smooth()` using formula 'y ~ x'

PlotScatter(tib, "distance_to_telomere", "RPE_0h_Ki67", xlimits = c(0, 50), 
            color_by = tib$RPE_wt_H3K9me3, smooth_line = T)
## `geom_smooth()` using formula 'y ~ x'

PlotScatter(tib, "distance_to_telomere", "RPE_0h_Ki67", xlimits = c(0, 50), 
            color_by = tib$RPE_wt_H3K27me3, smooth_seqnames = T)
## `geom_smooth()` using formula 'y ~ x'

# Repeat for t=1h
PlotScatter(tib, "distance_to_centromere", "RPE_1h_Ki67", xlimits = c(0, 50))

PlotScatter(tib, "distance_to_centromere", "RPE_1h_Ki67", xlimits = c(0, 50), 
            color_by = tib$RPE_wt_H3K27me3)

PlotScatter(tib, "distance_to_centromere", "RPE_1h_Ki67", xlimits = c(0, 50), 
            color_by = tib$RPE_wt_H3K9me3)

PlotScatter(tib, "distance_to_telomere", "RPE_1h_Ki67", xlimits = c(0, 50))

PlotScatter(tib, "distance_to_telomere", "RPE_1h_Ki67", xlimits = c(0, 50), 
            color_by = tib$RPE_wt_H3K27me3)

PlotScatter(tib, "distance_to_telomere", "RPE_1h_Ki67", xlimits = c(0, 50), 
            color_by = tib$RPE_wt_H3K9me3)

PlotScatter(tib, "distance_to_telomere", "RPE_1h_Ki67", xlimits = c(0, 50), 
            smooth_line = T)
## `geom_smooth()` using formula 'y ~ x'

PlotScatter(tib, "distance_to_telomere", "RPE_1h_Ki67", xlimits = c(0, 50), 
            color_by = tib$RPE_wt_H3K27me3, smooth_line = T)
## `geom_smooth()` using formula 'y ~ x'

PlotScatter(tib, "distance_to_telomere", "RPE_1h_Ki67", xlimits = c(0, 50), 
            color_by = tib$RPE_wt_H3K9me3, smooth_line = T)
## `geom_smooth()` using formula 'y ~ x'

PlotScatter(tib, "distance_to_telomere", "RPE_1h_Ki67", xlimits = c(0, 50), 
            color_by = tib$RPE_wt_H3K27me3, smooth_seqnames = T)
## `geom_smooth()` using formula 'y ~ x'

PlotScatter(tib, "distance_to_telomere", "RPE_0h_Ki67", xlimits = c(0, 50), 
            color_by = tib$RPE_wt_H3K27me3, smooth_line = T, facet_seqnames = T)
## `geom_smooth()` using formula 'y ~ x'

# t=0h
PlotScatterBinned(tib, "distance_to_telomere", "RPE_0h_Ki67", xlimits = c(0, 50),
                  smooth_line = T, limits_count = c(0, 100))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `geom_smooth()` using formula 'y ~ x'

PlotScatterBinned(tib, "distance_to_telomere", "RPE_0h_Ki67", xlimits = c(0, 50), 
                  color_by = tib$RPE_wt_H3K27me3, ylimits_col = c(-1.8, 1.8))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, "distance_to_telomere", "RPE_0h_Ki67", xlimits = c(0, 50), 
                  color_by = tib$RPE_wt_H3K9me3, ylimits_col = c(-2, 2))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

# t=1h
PlotScatterBinned(tib, "distance_to_telomere", "RPE_1h_Ki67", xlimits = c(0, 50),
                  smooth_line = T, limits_count = c(0, 100))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `geom_smooth()` using formula 'y ~ x'

PlotScatterBinned(tib, "distance_to_telomere", "RPE_1h_Ki67", xlimits = c(0, 50), 
                  color_by = tib$RPE_wt_H3K27me3, ylimits_col = c(-1.8, 1.8))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, "distance_to_telomere", "RPE_1h_Ki67", xlimits = c(0, 50), 
                  color_by = tib$RPE_wt_H3K9me3, ylimits_col = c(-2, 2))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

# Prepare smooth line of all telomere distances separately
tib %>%
  dplyr::select(seqnames, distance_to_telomere, RPE_0h_Ki67) %>%
  drop_na() %>%
  filter(distance_to_telomere < 55) %>%
  mutate(seqnames = factor(seqnames, chrom_sizes$seqnames)) %>%
  ggplot(aes(x = distance_to_telomere, y = RPE_0h_Ki67, col = seqnames)) +
  geom_smooth(method = "loess", se = F, span = 10) +
  geom_vline(xintercept = 0) +
  #geom_vline(xintercept = 25, linetype = "dashed") +
  geom_hline(yintercept = 0) +
  coord_cartesian(xlim = c(0, 50),
                  ylim = c(-0.5, 2)) +
  theme_bw() +
  theme(aspect.ratio = 1)  
## `geom_smooth()` using formula 'y ~ x'

6. Comparison with lamina changes during the cell cycle

# Load dynamic LADs
dynamic_lads <- read_rds("~/mydata/proj/tests/results/ts190509_RPE_shakeOff/ts190802_differential_analysis_BinsandLADs_RPE/LADs.rds")
dynamic_results <- read_rds("~/mydata/proj/tests/results/ts190509_RPE_shakeOff/ts190802_differential_analysis_BinsandLADs_RPE/LADs_results.rds")

dynamic_lads$result <- as.vector(dynamic_results[, 2])
## Loading required package: limma
## 
## Attaching package: 'limma'
## The following object is masked from 'package:BiocGenerics':
## 
##     plotMA
# Calculate change of Ki67 on dynamic LADs
ovl <- findOverlaps(gr_padamid_combined, dynamic_lads, select = "arbitrary")

tib_lads <- tib %>%
  dplyr::select(-start, -end, -contains("distance")) %>%
  mutate(match = findOverlaps(as(tib, "GRanges"), 
                              dynamic_lads, select = "arbitrary"),
         lad_result = dynamic_lads$result[match],
         lad_result = case_when(lad_result == 1 ~ "up",
                                lad_result == -1 ~ "down",
                                T ~ "stable"),
         lad_result = factor(lad_result, c("down", "stable", "up"))) %>%
  drop_na(match) %>%
  gather(key, value, contains("RPE")) %>%
  group_by(seqnames, key, match, lad_result) %>%
  dplyr::summarise(mean = mean(value, na.rm = T)) %>%
  ungroup() %>%
  spread(key, mean) %>%
  drop_na()
## `summarise()` has grouped output by 'seqnames', 'key', 'match'. You can override using the `.groups` argument.
# 1) Plot scores at early time points
tib_lads %>%
  dplyr::select(lad_result, RPE_0h_Ki67, RPE_1h_Ki67, RPE_3h_Ki67) %>%
  gather(key, value, -lad_result) %>%
  mutate(key = factor(key, unique(key))) %>%
  ggplot(aes(x = lad_result, y = value, col = lad_result)) + 
  geom_quasirandom() +
  #geom_boxplot(outlier.shape = NA, fill = NA, col = "black") +
  stat_summary(fun.data = quantiles, geom = "boxplot", fill = NA, col = "black") +
  geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
  facet_grid(. ~ key) +
  scale_color_brewer(palette = "Dark2", guide = F) +
  xlab("LAD class") +
  ylab("Ki67 t=0h") +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, hjust = 1))

# 2) Plot differences at dynamic lads
tib_lads %>%
  mutate(RPE_3h_Ki67 = RPE_3h_Ki67 - RPE_1h_Ki67,
         RPE_6h_Ki67 = RPE_6h_Ki67 - RPE_1h_Ki67,
         RPE_10h_Ki67 = RPE_10h_Ki67 - RPE_1h_Ki67) %>%
  dplyr::select(lad_result, RPE_3h_Ki67, RPE_6h_Ki67, RPE_10h_Ki67) %>%
  gather(key, value, -lad_result) %>%
  mutate(key = factor(key, unique(key))) %>%
  ggplot(aes(x = lad_result, y = value, col = lad_result)) + 
  geom_quasirandom() +
  #geom_boxplot(outlier.shape = NA, fill = NA, col = "black") +
  stat_summary(fun.data = quantiles, geom = "boxplot", fill = NA, col = "black") +
  geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
  facet_grid(. ~ key) +
  scale_color_brewer(palette = "Dark2", guide = F) +
  xlab("LAD class") +
  ylab("Ki67 difference with t=1h") +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, hjust = 1))

There is no simple correlation of increasing LaminB1 and decreasing Ki67 signals, and vice-versa, during the cell cycle. I think that this is due to the combination of aspects that determine the interaction patterns. In the other analyses, I can consistently find evidence that LaminB1 and Ki67 have opposite behaviour.

Conclusions

Several conclusions:

  • Ki67 interactions are widely spread early in G1 and mature towards small chromosomes.
  • Ki67 is enriched at the chromosome ends during mitosis / metaphase. Not sure why.
  • It seems that this enrichment scales with chromosome size. It could thus be related to the arms that escape the chromosomal mass, but also due to encapsulation of these ends.
  • While the correlation with H3K9me3-marked chromatin is lost, the anti correlation with H3K27me3 remains.
  • The enrichment on chromosomal arms can of course also be a technical artifact from the pA-DamID procedure.

Session info

sessionInfo()
## R version 4.0.5 (2021-03-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] limma_3.46.0         knitr_1.32           ggbeeswarm_0.6.0    
##  [4] caTools_1.18.2       corrr_0.4.3          GGally_2.1.1        
##  [7] RColorBrewer_1.1-2   rtracklayer_1.50.0   GenomicRanges_1.42.0
## [10] GenomeInfoDb_1.26.7  IRanges_2.24.1       S4Vectors_0.28.1    
## [13] BiocGenerics_0.36.1  forcats_0.5.1        stringr_1.4.0       
## [16] dplyr_1.0.5          purrr_0.3.4          readr_1.4.0         
## [19] tidyr_1.1.3          tibble_3.1.1         ggplot2_3.3.3       
## [22] tidyverse_1.3.1     
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-152                bitops_1.0-6               
##  [3] matrixStats_0.58.0          fs_1.5.0                   
##  [5] lubridate_1.7.10            httr_1.4.2                 
##  [7] tools_4.0.5                 backports_1.2.1            
##  [9] bslib_0.2.4                 utf8_1.2.1                 
## [11] R6_2.5.0                    vipor_0.4.5                
## [13] mgcv_1.8-35                 DBI_1.1.1                  
## [15] colorspace_2.0-0            withr_2.4.2                
## [17] tidyselect_1.1.0            compiler_4.0.5             
## [19] cli_2.4.0                   rvest_1.0.0                
## [21] Biobase_2.50.0              xml2_1.3.2                 
## [23] DelayedArray_0.16.3         labeling_0.4.2             
## [25] sass_0.3.1                  scales_1.1.1               
## [27] digest_0.6.27               Rsamtools_2.6.0            
## [29] rmarkdown_2.7               XVector_0.30.0             
## [31] pkgconfig_2.0.3             htmltools_0.5.1.1          
## [33] MatrixGenerics_1.2.1        highr_0.9                  
## [35] dbplyr_2.1.1                rlang_0.4.10               
## [37] readxl_1.3.1                rstudioapi_0.13            
## [39] farver_2.1.0                jquerylib_0.1.3            
## [41] generics_0.1.0              jsonlite_1.7.2             
## [43] BiocParallel_1.24.1         RCurl_1.98-1.3             
## [45] magrittr_2.0.1              GenomeInfoDbData_1.2.4     
## [47] Matrix_1.3-2                Rcpp_1.0.6                 
## [49] munsell_0.5.0               fansi_0.4.2                
## [51] lifecycle_1.0.0             stringi_1.5.3              
## [53] yaml_2.2.1                  SummarizedExperiment_1.20.0
## [55] zlibbioc_1.36.0             plyr_1.8.6                 
## [57] grid_4.0.5                  crayon_1.4.1               
## [59] lattice_0.20-41             splines_4.0.5              
## [61] Biostrings_2.58.0           haven_2.4.0                
## [63] hms_1.0.0                   pillar_1.6.0               
## [65] codetools_0.2-18            reprex_2.0.0               
## [67] XML_3.99-0.6                glue_1.4.2                 
## [69] evaluate_0.14               modelr_0.1.8               
## [71] vctrs_0.3.7                 cellranger_1.1.0           
## [73] gtable_0.3.0                reshape_0.8.8              
## [75] assertthat_0.2.1            xfun_0.22                  
## [77] broom_0.7.6                 beeswarm_0.3.1             
## [79] GenomicAlignments_1.26.0    ellipsis_0.3.1